#define _VARIADIC_MAX 10


#include <cmath>
#include <iterator>
#include <list>
#include <queue>
#include <omp.h>
#include <fstream>
#include "config.h"
#include "progressive_photon_map_renderer.h"
#include "photon_map.h"

#include "image_film.h"
#include "camera.h"
#include "random_number_generator_mt19937.h"
#include "permuted_halton.h"
#include "scene_data.h"
#include "volume_integrator.h"
#include "opencl_intersect_device.h"
#include "console_progress_reporter.h"

#include <ctime> 
#include "parallel.h"

PPMRenderer::PPMRenderer(Camera* c,ImageFilm* im,Sampler* s,photon_map_t* photon_map)
:camera_(c),image_(im),sampler_(s),photon_map_(photon_map)
{
	progress_reporter_ = new ConsoleProgressReporter();
	volume_integrator_ = new NullVolumeIntegrator();
}
PPMRenderer::~PPMRenderer()
{
	delete camera_;
	delete sampler_;
	delete image_;
	delete progress_reporter_;
	delete volume_integrator_;
	delete photon_map_;
}


const size_t ray_buffer_size = 1024 * 16;
const size_t ray_buffer_max_size = 1024 * 16;
const size_t max_queue_size = 8;

const size_t photon_buffer_size = 1024 * 16;
const size_t photon_buffer_max_size = 1024 * 16;

const int max_ray_depth = 6;

struct eye_path_t{
	camera_sample_t camera_sample;
	spectrum_t passthrough;
	//ray_differential_t ray;//size_t rayid,ray_t rx,ry;
	ray_t main_ray,rx,ry;
	int has_differential;
	int depth;
};

void PPMRenderer::RenderIteration(unsigned iteration)
{
	sampler_->ResetSamplePosition();
}

struct EstimateRadiance{
public:
	const photon_map_t *photon_map_;
	std::shared_ptr<Sampler> sampler_;
	std::shared_ptr<OpenCLIntersectDevice> intersect_device_;
	const cl_scene_info_t& cl_scene_info_;
	const Camera* camera_;
	const int task_id_ ;
	const VolumeIntegrator* volume_integrator_;

	mutable tbb::mutex& mutex_;
	mutable ImageFilm* film_;

	EstimateRadiance(std::shared_ptr<Sampler> sampler, const photon_map_t* photon_map,std::shared_ptr<OpenCLIntersectDevice> isect_device,
		const cl_scene_info_t &cl_scene_info,
		const Camera* camera,
		const VolumeIntegrator* volume_integrator,
		int task_id, ImageFilm* film,  
		tbb::mutex& mtx):
	sampler_(sampler),photon_map_(photon_map),intersect_device_(isect_device),
		cl_scene_info_(cl_scene_info), camera_(camera),  task_id_(task_id),film_(film),mutex_(mtx),volume_integrator_(volume_integrator)
	{
	}
	void operator()()const
	{
		bool has_more_sample = true;
		int ray_traced = 0;

		std::vector<std::pair<camera_sample_t,spectrum_t> > result;

		
		std::list<eye_path_t> todo_eye_paths;
		std::queue<std::list<eye_path_t>> todo_eye_paths_queue;
		std::shared_ptr<std::vector<ray_t> > ray_buffer = intersect_device_->CreateRayBuffer(ray_buffer_size);

		ray_buffer->clear();
		ray_buffer->reserve(ray_buffer_size);

		std::list<eye_path_t>::iterator todo_eye_path_iterator = todo_eye_paths.begin();
		//lphoton data
		close_photon_t close_photon_data_store[MAX_CLOSE_PHOTON_LOOKUP];
		//final gather data
		close_photon_t photon_buffer[n_indir_sample_photons];
		vector3f_t photon_dirs [n_indir_sample_photons];
		sampler_->ResetSamplePosition();
		Seed seed;
		init_rng((31 * (1+photon_map_->progressive_iteration) ) ^ (task_id_ * 29),&seed);
		while(has_more_sample || (!todo_eye_paths_queue.empty() || !todo_eye_paths.empty()))//do eye pass
		{
			//
			camera_sample_t sample;
			ray_differential_t ray;
			float ray_weight;
			
			todo_eye_path_iterator = todo_eye_paths.begin();
			while( todo_eye_path_iterator != todo_eye_paths.end() )
			{
				if (todo_eye_path_iterator->depth < max_ray_depth)
				{
					todo_eye_path_iterator->depth ++;
					ray_buffer->push_back(todo_eye_path_iterator->main_ray);
					todo_eye_path_iterator++;
				}
				else
				{
					todo_eye_path_iterator = todo_eye_paths.erase(todo_eye_path_iterator);
				}
			}
			while(ray_buffer->size() < ray_buffer_size && has_more_sample)
			{
				if(sampler_->GetNextSample(&sample))
				{
					camera_->GenerateRay(sample, &ray, &ray_weight);
					ray_t r;
					rassign(r,ray);
					ray_buffer->push_back(r);
					eye_path_t todo_eye_path;
					todo_eye_path.main_ray = r;
					todo_eye_path.rx = ray.rx;
					todo_eye_path.ry = ray.rx;
					todo_eye_path.has_differential = ray.has_differential;
					todo_eye_path.camera_sample = sample;
					todo_eye_path.depth = 0;
					todo_eye_path.passthrough = spectrum_t(1,1,1);
					todo_eye_paths.push_back(todo_eye_path);
				}
				else
					has_more_sample = false;
			}
			std::vector<intersection_t> intersections;
			if(ray_buffer->size() > 0)
			{
				todo_eye_paths_queue.push(std::list<eye_path_t>());
				todo_eye_paths_queue.back().swap(todo_eye_paths);
				intersect_device_->PushRayBuffer(ray_buffer);
			}
		
			if(ray_buffer->size() == 0)
			{
				int a_breeak = 0; 
			}
			if(intersect_device_->PopIntersectionBuffer(&ray_buffer, &intersections ))
			{
				todo_eye_paths_queue.front().swap(todo_eye_paths);
				todo_eye_paths_queue.pop();
			}
			todo_eye_path_iterator = todo_eye_paths.begin();
			for (size_t i = 0;i < intersections.size() && !todo_eye_paths.empty(); ++i)
			{
				(*ray_buffer)[i].maxt = sqrt(distance_squared(intersections[i].dg.p,(*ray_buffer)[i].o)) - (*ray_buffer)[i].mint;
				spectrum_t transmittance;
				spectrum_t volume_li = volume_integrator_->Li((*ray_buffer)[i],&transmittance);
				if(intersections[i].primitive_idx == 0xffffffff) //miss
				{
					vmul(volume_li,volume_li,todo_eye_path_iterator->passthrough);
					if(!color_is_black(volume_li))
					{
						result.push_back(std::make_pair(todo_eye_path_iterator->camera_sample,volume_li));
					}
					todo_eye_path_iterator = todo_eye_paths.erase(todo_eye_path_iterator);
				}
				else //hit
				{
					/*
					{
						ray_t test;
						point3f_t o={278.00000 ,273.00000 ,-800.00000 };
						vector3f_t d={-0.18108650 ,-0.012522597 ,0.98338735};
						intersection_t t;
						rinit(test,o,d);
						intersect(&cl_scene_info_.accelerator_data[0],&cl_scene_info_.shape_data[0],&cl_scene_info_.primitives[0],
							cl_scene_info_.primitive_count,&test,&t);
					}
					*/
					//cur_ray == (*ray_buffer)[i];
					bsdf_t bsdf;
					ray_differential_t cur_ray, out_ray;
					spectrum_t passthrough;
					rassign(cur_ray,todo_eye_path_iterator->main_ray);
					rassign(cur_ray.rx, todo_eye_path_iterator->rx);
					rassign(cur_ray.ry, todo_eye_path_iterator->ry);
					cur_ray.has_differential = todo_eye_path_iterator->has_differential;
					intersection_get_bsdf(&intersections[i],cl_scene_info_,
						&cur_ray,&bsdf);

					;
					
					bool specular_bounce = false;
					specular_reflect(&cur_ray,&bsdf,&seed,&intersections[i],cl_scene_info_,&passthrough,
                            &out_ray);
					if (!color_is_black(passthrough))
					{
						eye_path_t eye_path = *todo_eye_path_iterator;
						vmul(eye_path.passthrough,eye_path.passthrough,passthrough);
						vmul(eye_path.passthrough,eye_path.passthrough,transmittance);
						vadd(eye_path.passthrough,eye_path.passthrough,volume_li);
						rassign(eye_path.main_ray , out_ray);
						rassign(eye_path.rx, out_ray.rx);
						rassign(eye_path.ry, out_ray.ry);
						eye_path.has_differential = out_ray.has_differential;
						eye_path.main_ray.mint = intersections[i].ray_epsilon;
						todo_eye_paths.push_back(eye_path);
						specular_bounce = true;
					}
					
					specular_transmit(&cur_ray,&bsdf,&seed,&intersections[i],cl_scene_info_,&passthrough,
                                &out_ray);
					if (!color_is_black(passthrough))
					{
						eye_path_t eye_path = *todo_eye_path_iterator;
						vmul(eye_path.passthrough,eye_path.passthrough,passthrough);
						vmul(eye_path.passthrough,eye_path.passthrough,transmittance);
						vadd(eye_path.passthrough,eye_path.passthrough,volume_li);

						rassign(eye_path.main_ray , out_ray);
						rassign(eye_path.rx, out_ray.rx);
						rassign(eye_path.ry, out_ray.ry);
						eye_path.has_differential = out_ray.has_differential;
						eye_path.main_ray.mint = intersections[i].ray_epsilon;
						todo_eye_paths.push_back(eye_path);
						specular_bounce = true;
					}
					if(specular_bounce)
					{
						todo_eye_path_iterator = todo_eye_paths.erase(todo_eye_path_iterator);
					}
					else
					{
						if (int(todo_eye_path_iterator->camera_sample.image_x) == 135 &&
							int(todo_eye_path_iterator->camera_sample.image_y) == 507)
						{
							int debug_break = 0;
						}
						vector3f_t wo;
						vneg(wo,cur_ray.d);
						//compute emitted light if ray hit an area light
						spectrum_t l;
						spectrum_t li_val;
						vclr(li_val);
						vclr(l);
						intersection_le(&intersections[i],cl_scene_info_,&wo,&l);
						vadd(li_val,li_val,l);
						point3f_t p = bsdf.dg_shading.p;
						normal3f_t n = bsdf.dg_shading.nn;
						uniform_sample_all_lights(cl_scene_info_,
								&seed,&p,&n,&wo,&bsdf,&l);
						vadd(li_val,li_val,l);

						//#define DIRECT_LIGHTING
#ifndef DIRECT_LIGHTING
						photon_map_lphoton(photon_map_,
							photon_map_->caustic_map_lookup,photon_map_->n_caustic_paths,
							photon_map_->n_lookup,&bsdf,&intersections[i],&wo,photon_map_->max_dist_squared
							,cl_scene_info_,&seed,&l,close_photon_data_store);
						vadd(li_val,li_val,l);
						if (photon_map_->final_gather)
						{
#if 1
							photon_map_final_gather(photon_map_,
								cl_scene_info_,&seed,p,n,wo,&bsdf,&l,
								photon_buffer,
								photon_dirs);
#else
							normal3f_t nn;
							vassign(nn,n);
							if(vdot(nn,cur_ray.d) > 0.f) vneg(nn,n);
							radiance_photon_process_data_t rpd;
							radiance_photon_process_data_init(&rpd,&p,&nn);
							float md2 = FLT_MAX;
							kd_tree_lookup( photon_map->radiance_map,p,&rpd, radiance_photon_process,md2);
							if(rpd.photon)
								vadd(li_val,li_val,rpd.photon->lo);
#endif
						}
						else
							photon_map_lphoton(photon_map_,
							photon_map_->indirect_map_lookup,
							photon_map_->n_indirect_paths,
							photon_map_->n_lookup,&bsdf,&intersections[i],&wo,photon_map_->max_dist_squared
							,cl_scene_info_,&seed,&l,close_photon_data_store);
#endif
						vadd(li_val,li_val,l);
						vmul(li_val,li_val,todo_eye_path_iterator->passthrough);
						spectrum_t color = li_val;
						vmul(color,color,transmittance);
						vadd(color,color,volume_li);
						result.push_back(std::make_pair(todo_eye_path_iterator->camera_sample,color));
						todo_eye_path_iterator = todo_eye_paths.erase(todo_eye_path_iterator);
					}
				}
			{
				tbb::mutex::scoped_lock lock(mutex_);
				for(auto it = result.begin();it != result.end(); ++it)
				{
					film_->AddSample((*it).first,(*it).second);
				}
			}
				result.clear();

			}

			ray_buffer->clear();
		}

	}
};
void PPMRenderer::Render(const scene_info_memory_t& scene_info_mem)
{
	RandomNumberGeneratorMT19937 *rng = new RandomNumberGeneratorMT19937(rand() << 16 | rand());

	scene_info_memory_t scene_info = scene_info_mem;
	cl_scene_info_t cl_scene_info = as_cl_scene_info(scene_info);
	OpenCLIntersectDevice::InitializeKernel(scene_info);
	int iteration = 0;//start_iteration

	int n = iteration + 1;
	float ppm_ri_r1 = 1.f / n;
	for(int i = 1;i < n;++i)
	{
		ppm_ri_r1 *= (i+photon_map_->alpha)/(i);
	}
	photon_map_->max_dist_squared *= ppm_ri_r1;
	while(true)
	{
		clock_t t0;
		t0 = clock();

		InitializePhotonMap(photon_map_,scene_info,rng);
		std::list<std::shared_ptr<EstimateRadiance> > all_estimators;

		size_t total_sample_count = sampler_->GetTotalSampleCount();

		tbb::task_group all_task_group;
		int total_task_count = 8*get_default_num_threads();
		tbb::mutex mutex;
		for(int i = 0;i < total_task_count;++i)
		{
			std::shared_ptr<Sampler> sub_sampler = sampler_->GetSubSampler(i,total_task_count);
			std::shared_ptr<OpenCLIntersectDevice> intersect_device
				(new OpenCLIntersectDevice(scene_info,ray_buffer_size, max_queue_size));
			if (sub_sampler)
			{
				all_estimators.push_back(
					std::make_shared<EstimateRadiance>(
					sub_sampler,photon_map_,intersect_device,cl_scene_info,camera_,volume_integrator_,
					i,
					image_,mutex
					));
				all_task_group.run(tbb::make_task(*all_estimators.back()));
			}
		}

		all_task_group.wait();

		photon_map_destroy(photon_map_);
		image_->WriteImage(iteration);
		image_->ClearSamples();
		
		
		iteration ++;
		photon_map_->max_dist_squared *= ((iteration + photon_map_->alpha)/(iteration + 1));
		photon_map_->progressive_iteration = iteration;
		printf("iteration:%d, max_dist_squared:%.3f\n",iteration,photon_map_->max_dist_squared);
	}
	delete rng;
}

struct photon_path_t{
	spectrum_t alpha;
	ray_t ray;
	unsigned int depth;
	bool is_specular;
	bool is_caustic_deposited;
	bool is_indirect_deposited;
	bool is_direct_deposited;
};

struct BuildPhotonMapResult{
	std::vector<photon_t> caustic_photons;
	std::vector<photon_t> indirect_photons;
	std::vector<photon_t> direct_photons;
	std::vector<radiance_photon_t> radiance_photons;

	std::vector<spectrum_t> rp_reflectances,rp_transmittances;
	int n_shot;
	size_t caustic_done_paths_;
	size_t indirect_done_paths_;
	size_t direct_done_paths_;

	BuildPhotonMapResult():n_shot(0),caustic_done_paths_(0),indirect_done_paths_(0),direct_done_paths_(0){}
};
struct BuildPhotonMap{
	const photon_map_t *photon_map_;
	const VolumeIntegrator* volume_integrator_;
	const cl_scene_info_t& cl_scene_info_;
	std::shared_ptr<OpenCLIntersectDevice> intersect_device_;
	const int task_id_;
	const int start_total_shot_;
	const size_t n_caustic_photons_;
	const size_t n_indirect_photons_;

	static const int halton_dimension = 8;
	mutable std::shared_ptr<BuildPhotonMapResult>  result_photons_;

	tbb::mutex &mutex_;

	BuildPhotonMap(
		const photon_map_t *photon_map,
		const cl_scene_info_t& cl_scene_info,
		const VolumeIntegrator* volume_integrator,
		std::shared_ptr<OpenCLIntersectDevice> intersect_device,
		const int task_id,const int start_total_shot,const size_t n_caustic_photons,const size_t n_indirect_photons,
		std::shared_ptr<BuildPhotonMapResult> result_photons,
		tbb::mutex& mutex
		):photon_map_(photon_map),cl_scene_info_(cl_scene_info),intersect_device_(intersect_device),task_id_(task_id),
		start_total_shot_(start_total_shot),n_caustic_photons_(n_caustic_photons),n_indirect_photons_(n_indirect_photons),
		result_photons_(result_photons),mutex_(mutex),volume_integrator_(volume_integrator)
	{
	}
	void operator()()const
	{
		using namespace std;

		//thread local photons
		std::vector<photon_t> caustic_photons;
		std::vector<photon_t> indirect_photons;
		std::vector<photon_t> direct_photons;
		std::vector<radiance_photon_t> radiance_photons;
		std::vector<spectrum_t> rp_reflectances,rp_transmittances;
		size_t caustic_done_paths = 0;
		size_t indirect_done_paths = 0;
		size_t direct_done_paths = 0;
	
		caustic_photons.reserve(n_caustic_photons_);
		indirect_photons.reserve(n_indirect_photons_);

		int n_shot=0;
		bool caustics_done = (n_caustic_photons_ == 0);
		bool indirect_done = (n_indirect_photons_ == 0);
		//compute light power cdf for photon shooting
		int n_lights = int(cl_scene_info_.lght_count);
		float *lights_power = (float*)alloca	(n_lights * sizeof(float));
		float *light_cdf = (float*)alloca((n_lights+1) * sizeof(float));
		spectrum_t lpower;
		for (int i = 0;i < n_lights; ++i)
		{
			light_power(cl_scene_info_.lghts+i,cl_scene_info_,&lpower);
			lights_power[i] = spectrum_y(&lpower);
		}
		float total_power;
		compute_step_1d_cdf(lights_power,n_lights,&total_power,light_cdf);
		//radiance photon reflectance arrays


		Seed seed;
		init_rng(37 * (1+photon_map_->progressive_iteration) ^ (43 * task_id_),&seed);
		RandomNumberGeneratorMT19937 rng((53*(1+photon_map_->progressive_iteration)) ^ (101 * task_id_));
		permuted_halton16_t halton ;
		permuted_halton16_initialize(&halton,rng);

		//static const int buffer_size = 1 * 16 * 128;
		std::shared_ptr<std::vector<ray_t> > photon_rays = intersect_device_->CreateRayBuffer(photon_buffer_size);
		photon_rays->clear();
		std::list<photon_path_t> todo_photon_paths;
		std::queue<std::list<photon_path_t>> todo_photon_paths_queue;
		std::queue<std::list<photon_path_t>> resident_photon_paths_queue;
		std::list<photon_path_t>::iterator todo_photon_path_iterator = todo_photon_paths.begin();

		ray_differential_t the_photon_ray_differential;

		while (!caustics_done || !indirect_done )//todo: add not photon_rays.empty
		{
			if(!resident_photon_paths_queue.empty())
			{
				todo_photon_paths.swap(resident_photon_paths_queue.front());
				resident_photon_paths_queue.pop();
			}
			todo_photon_path_iterator = todo_photon_paths.begin();
			assert(photon_rays->size() == 0 && todo_photon_paths.size() <= photon_buffer_size);
			while(photon_rays->size() < photon_buffer_size && todo_photon_path_iterator != todo_photon_paths.end())
			{
				if(todo_photon_path_iterator->depth < photon_map_->max_specular_depth)
				{
					photon_rays->push_back(todo_photon_path_iterator->ray);
					todo_photon_path_iterator++;
				}
				else
				{
					todo_photon_path_iterator = todo_photon_paths.erase(todo_photon_path_iterator);
				}
			}
			while(photon_rays->size() < photon_buffer_size)
			{
				if(resident_photon_paths_queue.empty() && (caustic_photons.size() < n_caustic_photons_ 
				|| indirect_photons.size() < n_indirect_photons_))
				{
					n_shot ++;
					float u[5];

					ray_t photon_ray;
					spectrum_t alpha;

					halton16_sample(&halton,start_total_shot_ + n_shot, u , 5);
					//choose light of shoot photon from
					float lpdf;

					int lnum = (int)floor(sample_step_1d(lights_power,light_cdf,
						total_power,n_lights,u[0],&lpdf));
					lnum = min(lnum, (int)n_lights - 1);

					GLOBAL const light_info_t* light = cl_scene_info_.lghts +lnum;
					float pdf;

					normal3f_t nl;
					light_ray_sample_l(light,cl_scene_info_,u[1 ],u[2 ],u[3],u[4],
						&photon_ray,&nl,&pdf,&alpha);
					if (pdf == 0.f || color_is_black(alpha))
					{
						continue;
					}
					photon_path_t photon_path;
					vsmul(photon_path.alpha,(fabs(vdot(nl,photon_ray.d))/(pdf*lpdf)),alpha);
					photon_path.depth = 0;
					photon_path.ray = photon_ray;
					photon_path.is_specular = true;
					photon_path.is_caustic_deposited = false;
					photon_path.is_indirect_deposited = false;
					photon_path.is_direct_deposited = false;
					photon_rays->push_back(photon_ray);
					todo_photon_paths.push_back(photon_path);
				}
				else
				{
					if(!resident_photon_paths_queue.empty())
					{
						std::list<photon_path_t> &resident_photon_paths = resident_photon_paths_queue.front();
						auto resident_photon_path_iterator = resident_photon_paths.begin();
						while(photon_rays->size() < photon_buffer_size && resident_photon_path_iterator != resident_photon_paths.end())
						{
							if(resident_photon_path_iterator->depth < photon_map_->max_specular_depth)
							{
								photon_rays->push_back(resident_photon_path_iterator->ray);
								todo_photon_paths.push_back(*resident_photon_path_iterator);
							}
							resident_photon_path_iterator = resident_photon_paths.erase(resident_photon_path_iterator);
						}
						if(resident_photon_path_iterator == resident_photon_paths.end())
						{
							resident_photon_paths_queue.pop();
						}
					}		
					if(resident_photon_paths_queue.empty()) break;
				}
			}
			std::vector<intersection_t> intersections;
			if(!todo_photon_paths.empty())
			{
				todo_photon_paths_queue.push(std::list<photon_path_t>());
				todo_photon_paths_queue.back().swap(todo_photon_paths);
				intersect_device_->PushRayBuffer(photon_rays);
			}

			do{
				if(intersect_device_->PopIntersectionBuffer(&photon_rays,&intersections))
				{
					todo_photon_paths_queue.front().swap(todo_photon_paths);
					todo_photon_paths_queue.pop();
				}
				else
				{
					todo_photon_paths.clear();
				}
				assert(intersections.size() == todo_photon_paths.size());
				todo_photon_path_iterator = todo_photon_paths.begin();
				for(size_t i = 0;i < intersections.size(); ++i)
				{

					(*photon_rays)[i].maxt = sqrt(distance_squared(intersections[i].dg.p,(*photon_rays)[i].o)) - (*photon_rays)[i].mint;
					if(intersections[i].primitive_idx == 0xffffffff)
					{
						todo_photon_path_iterator = todo_photon_paths.erase(todo_photon_path_iterator);
					}
					else
					{
						vector3f_t wo ; 
						vneg(wo,todo_photon_path_iterator->ray.d);
						spectrum_t transmittance = volume_integrator_->Transmittance((*photon_rays)[i]);
						bsdf_t photon_bsdf;
						rassign(the_photon_ray_differential, todo_photon_path_iterator->ray)
							intersection_get_bsdf(&intersections[i],cl_scene_info_,&the_photon_ray_differential,
							&photon_bsdf);
						BxDFType specular_type = 
							(BxDFType)(BSDF_REFLECTION|BSDF_TRANSMISSION|BSDF_SPECULAR);
						bool has_non_specular =	photon_bsdf.n_bxdfs > bsdf_num_components(&photon_bsdf,specular_type);
						todo_photon_path_iterator->depth++;
						if (todo_photon_path_iterator->depth > 0 && has_non_specular)
						{
							//deposit photon at surface
							photon_t photon;
							spectrum_t color;
							vmul(color,transmittance,todo_photon_path_iterator->alpha);

							photon_init(&photon,&intersections[i].dg.p,&color,&wo);
							bool deposited = false;
							if(todo_photon_path_iterator->depth == 1)
							{
								if(!indirect_done &&  photon_map_->final_gather)
								{
									direct_photons.push_back(photon);
									deposited = true;

									if(!todo_photon_path_iterator->is_direct_deposited)
									{
										(direct_done_paths) ++;
										todo_photon_path_iterator->is_direct_deposited = true;
									}
								}
							}
							else
							{
								if(todo_photon_path_iterator->is_specular)
								{
									//if (!*caustics_done_)
									{
										caustic_photons.push_back(photon);
										deposited = true;
										
										if(!todo_photon_path_iterator->is_caustic_deposited)
										{
											(caustic_done_paths) ++;
											todo_photon_path_iterator->is_caustic_deposited = true;
										}
									}
								}
								else
								{
									//if (!*indirect_done_)
									{
										indirect_photons.push_back(photon);
										deposited = true;
										if(!todo_photon_path_iterator->is_indirect_deposited)
										{
											(indirect_done_paths) ++;
											todo_photon_path_iterator->is_indirect_deposited = true;
										}
									}
								}
							}
							if(deposited && photon_map_->final_gather && random_float(&seed) < photon_map_->rr_threshold)
							{
								//store data for radiance photon
								normal3f_t n = intersections[i].dg.nn;
								radiance_photon_t radiance_photon;
								spectrum_t rho_r,rho_t;
								if(vdot(n,todo_photon_path_iterator->ray.d) > 0.f)vneg(n,n);
								radiance_photon_init(&radiance_photon,&(intersections[i].dg.p),(&n));
								bsdf_rho_hh(&photon_bsdf,&seed,BSDF_ALL_REFLECTION,&rho_r);
								bsdf_rho_hh(&photon_bsdf,&seed,BSDF_ALL_TRANSMISSION,&rho_t);
								radiance_photons.push_back(radiance_photon);
								rp_reflectances.push_back(rho_r);
								rp_transmittances.push_back(rho_t);
							}
						}
						//

						//sample new photon ray direction
						vector3f_t wi;
						float pdf;
						BxDFType flags;

						//compute new photon weight and possibly terminate with rr
						float u1 = random_float(&seed);
						float u2 = random_float(&seed);
						float u3 = random_float(&seed);
						spectrum_t fr;
						bsdf_sample_f(&photon_bsdf,&wo,&wi,u1,u2,u3,&pdf,BSDF_ALL,&flags,&fr);
						if(color_is_black(fr) || pdf == 0.f)
						{
							todo_photon_path_iterator = todo_photon_paths.erase(todo_photon_path_iterator);
							continue;
						}
						todo_photon_path_iterator->is_specular = todo_photon_path_iterator->is_specular && ((flags & BSDF_SPECULAR) != 0);
						if(indirect_done && !todo_photon_path_iterator->is_specular)
						{
							todo_photon_path_iterator = todo_photon_paths.erase(todo_photon_path_iterator);
							continue;
						}

						spectrum_t anew;
						vmul(anew, todo_photon_path_iterator->alpha,fr);
						vsmul(anew,fabs(vdot(wi,photon_bsdf.dg_shading.nn))/pdf,anew);
						float continue_prob = min(1.f,spectrum_y(&anew)/spectrum_y(&todo_photon_path_iterator->alpha));
						if (random_float(&seed) > continue_prob  )
						{
							todo_photon_path_iterator = todo_photon_paths.erase(todo_photon_path_iterator);
							continue;
						}
						vsmul(todo_photon_path_iterator->alpha,1.f/continue_prob,anew);
						vmul(todo_photon_path_iterator->alpha,todo_photon_path_iterator->alpha,transmittance);

						rinit(todo_photon_path_iterator->ray,intersections[i].dg.p,wi);
						todo_photon_path_iterator->ray.mint = intersections[i].ray_epsilon;
						todo_photon_path_iterator++;
					}
				}
				if(!todo_photon_paths.empty())
				{
					resident_photon_paths_queue.push(std::list<photon_path_t>());
					resident_photon_paths_queue.back().swap(todo_photon_paths);
				}
			}
			while((caustic_photons.size() >= n_caustic_photons_ 
				&& indirect_photons.size() >= n_indirect_photons_)
				&& !todo_photon_paths_queue.empty());

			if (caustic_photons.size() >= n_caustic_photons_ && !caustics_done && resident_photon_paths_queue.empty())
			{
				//*caustic_done_paths_ = n_paths;//result_photons_->caustic_photons.size();
				caustics_done = true;
				//printf("%d:caustic_done:%d caustic_photon_size:%d\n",task_id_,n_shot,result_photons_->caustic_photons.size());
			}
			if (indirect_photons.size() >= n_indirect_photons_ && !indirect_done && resident_photon_paths_queue.empty())
			{
				//*indirect_done_paths_ = n_paths;//result_photons_->indirect_photons.size();
				indirect_done = true;
				//printf("%d:indirect_done:%d indirect_photon_size:%d\n",task_id_,n_shot,result_photons_->indirect_photons.size());
			}
			//give up if we are not storing enough photons
			if (n_shot > 5000000 && 
				(unsuccessful((unsigned)n_caustic_photons_,
				(unsigned)caustic_photons.size(),
				(unsigned)photon_buffer_size)) && 
				(unsuccessful((unsigned)n_indirect_photons_,
				(unsigned)indirect_photons.size(),
				(unsigned)photon_buffer_size)))
			{
				break;
			}
			photon_rays->clear();
		}

			{
				tbb::mutex::scoped_lock lock(mutex_);
				result_photons_->n_shot += n_shot;
				result_photons_->caustic_done_paths_ += caustic_done_paths;
				result_photons_->indirect_done_paths_ += indirect_done_paths;
				result_photons_->direct_done_paths_ += direct_done_paths;

				std::copy(indirect_photons.begin(),indirect_photons.end(),
					std::back_inserter(result_photons_->indirect_photons));
				indirect_photons.resize(0);
				std::copy(caustic_photons.begin(),caustic_photons.end(),
					std::back_inserter(result_photons_->caustic_photons));
				caustic_photons.resize(0);
				std::copy(radiance_photons.begin(),radiance_photons.end(),
					std::back_inserter(result_photons_->radiance_photons)
					);
				radiance_photons.resize(0);
				std::copy(direct_photons.begin(),direct_photons.end(),
					std::back_inserter(result_photons_->direct_photons)
					);
				direct_photons.resize(0);
				std::copy(rp_reflectances.begin(),rp_reflectances.end(),
					std::back_inserter(result_photons_->rp_reflectances)
					);
				rp_reflectances.resize(0);
				std::copy(rp_transmittances.begin(),rp_transmittances.end(),
					std::back_inserter(result_photons_->rp_transmittances)
					);
			}
	}

	static bool unsuccessful(int needed, int found, size_t shot) {
		return (found < needed &&
			(found == 0 || found < shot / 1024));
	}
};
void PPMRenderer::InitializePhotonMap(photon_map_t *photon_map, const scene_info_memory_t& scene_info_memory, 
	RandomNumberGeneratorMT19937 *rng)
{
	photon_map->indirect_map_lookup = NULL;
	photon_map->n_lookup = min(photon_map->n_lookup,MAX_CLOSE_PHOTON_LOOKUP);
	photon_map->max_specular_depth = min(photon_map->max_specular_depth,max_ray_depth);
	if (scene_info_memory.lghts.empty() ) return ;

	cl_scene_info_t cl_scene_info = as_cl_scene_info(scene_info_memory);
	clock_t t0 = clock();

	//
	using namespace std;
	
	tbb::task_group tg;
	tbb::mutex mutex;
	int task_count = 4*get_default_num_threads();
	std::list<std::shared_ptr<BuildPhotonMap> > all_build_photon_maps;
	std::shared_ptr<BuildPhotonMapResult> result(new BuildPhotonMapResult());
	for (int i = 0;i < task_count; ++i)
	{
		std::shared_ptr<OpenCLIntersectDevice> intersect_device
			(new OpenCLIntersectDevice(scene_info_memory,photon_buffer_size, max_queue_size));
		if (result && intersect_device)
		{
			all_build_photon_maps.push_back(std::make_shared<BuildPhotonMap>(photon_map_,cl_scene_info,volume_integrator_,
				intersect_device,
				i,photon_map_->total_photons + max(0,((photon_map_->n_caustic_photons + photon_map_->n_indirect_photons) * i)/task_count),
				max(0,photon_map_->n_caustic_photons/task_count),max(0,photon_map_->n_indirect_photons/task_count),
				result,mutex));
			tg.run(*all_build_photon_maps.back());
		}
	}
	tg.wait();

	bool indirect_map_inited = false;
	bool caustic_map_inited = false;

	std::shared_ptr<std::vector<photon_t> > caustic_photons(new std::vector<photon_t> ());
	caustic_photons->swap(result->caustic_photons);
	std::shared_ptr<std::vector<photon_t> > indirect_photons(new std::vector<photon_t> ());
	indirect_photons->swap(result->indirect_photons);

	if ( photon_map->n_caustic_photons > 0 && caustic_photons->size() >= photon_map->n_caustic_photons )
	{
		photon_map->n_caustic_paths =static_cast<unsigned int>(result->caustic_done_paths_);
		photon_map->caustic_map_lookup = new HybridLookupAccelerator<photon_t>
			(photon_map->n_lookup, photon_map->max_dist_squared,caustic_photons);
		caustic_map_inited = true;
	}
	if (photon_map->n_indirect_photons > 0 && indirect_photons->size() >= photon_map->n_indirect_photons )
	{
		photon_map->n_indirect_paths = static_cast<unsigned int>(result->indirect_done_paths_);
		//avoid incrementally search nearby photon in dead loop during search indirect photons
		if(photon_map->final_gather)
			photon_map->indirect_map_lookup = new HybridLookupAccelerator<photon_t>(indirect_photons);  
		else
			photon_map->indirect_map_lookup = 
			new HybridLookupAccelerator<photon_t>(photon_map->n_lookup, photon_map->max_dist_squared,indirect_photons);
		indirect_map_inited = true;
	}
	progress_reporter_->AddProgressInfo(
		photon_tracing_progress_reporter(caustic_photons->size(),indirect_photons->size(),result->direct_photons.size()));
	photon_map->total_photons += result->n_shot;
	if (photon_map->final_gather)
	{
		//precompute radiance at a subset of the photons
		std::shared_ptr<std::vector<photon_t> > direct_photons(new std::vector<photon_t>() );
		direct_photons->swap(result->direct_photons);
		HybridLookupAccelerator<photon_t>* direct_map = new 
			HybridLookupAccelerator<photon_t>(photon_map->n_lookup,photon_map->max_dist_squared,direct_photons);
			;
			t0 = clock();
		tbb::parallel_for(tbb::blocked_range<size_t>(0,result->radiance_photons.size()),
			[&](const tbb::blocked_range<size_t>& r)
			{
				for(size_t i = r.begin(); i != r.end(); ++i)
				{
					radiance_photon_t &rp = (result->radiance_photons)[i];
					const spectrum_t& rho_r = result->rp_reflectances[i];
					const spectrum_t& rho_t = result->rp_transmittances[i];
					spectrum_t E,E1,E2;
					point3f_t p = rp.p;
					normal3f_t n = rp.n;
					if (!color_is_black(rho_r))
					{
						vclr(E);
						vclr(E1);
						vclr(E2);
						estimate_e(photon_map,(direct_map),result->direct_done_paths_,p,n,&E);
						if(indirect_map_inited)
							estimate_e(photon_map,(photon_map->indirect_map_lookup),
							photon_map->n_indirect_paths,p,n,&E1);
						if(caustic_map_inited)
							estimate_e(photon_map,(photon_map->caustic_map_lookup),
							photon_map->n_caustic_paths,p,n,&E2);
						vadd(E,E,E1);
						vadd(E,E,E2);
						vmul(E,E,rho_r);
						vsmul(E,INV_PI,E);
						vadd(rp.lo,rp.lo,E);
					}
					if(!color_is_black(rho_t))
					{
						vclr(E);
						vclr(E1);
						vclr(E2);
						vneg(n,n);
						estimate_e(photon_map,(direct_map),result->direct_done_paths_,p,n,&E);
						if(indirect_map_inited)
							estimate_e(photon_map,(photon_map->indirect_map_lookup),
							photon_map->n_indirect_paths,p,n,&E1);
						if(caustic_map_inited)
							estimate_e(photon_map,(photon_map->caustic_map_lookup),
							photon_map->n_caustic_paths,p,n,&E2);
						vadd(E,E,E1);
						vadd(E,E,E2);
						vmul(E,E,rho_t);
						vsmul(E,INV_PI,E);
						vadd(rp.lo,rp.lo,E);
					}
				}
		});
		delete direct_map;
		std::shared_ptr<std::vector<radiance_photon_t> > radiance_photons(new std::vector<radiance_photon_t>);
		radiance_photons->swap(result->radiance_photons);
		photon_map->radiance_map_lookup = new HybridLookupAccelerator<radiance_photon_t>(photon_map->n_lookup,photon_map->max_dist_squared,radiance_photons);
	}
}



